Estudio Datos Vino. Depuración de variables

En este documento se exploran distintas posibilidades de estudio y depuración de un conjunto de datos con aplicación concreta al dataset de venta de vinos que se encuentra en la carpeta de Datos del material en formato cvs. Dado que la fase de “limpieza” de los datos es a crucial de cara a la obtención de modelos predictivos lo más libres posible de sesgos y problemas de estimación, es conveniente prestar atención e invertir el tiempo necesario en este proceso quizá explorando varias alternativas de las disponibles en base a la naturaleza de las distribuciones de las variables implicadas.

La fase de depuración de los datos no suele atender a “recetas” pues con cada dataset, el cuerpo nos pedirá cosas distintas. Sin embargo, trataremos de generar un proceso lo más automático posible con los aspectos generales que seguro tenemos que atacar para poder tenerlo como guía de actuación.

Estudio descriptivo de datos sobre venta de vinos

En primer lugar, importamos las dos librerías por excelencia, numpy y pandas para trabajar con conjuntos de datos en python, y leemos los datos de entrada consignando la ruta adecuada a nuestra carpeta (evidentemente hay que cambiar la ruta en cada caso, dudo que tengáis carpetar Guille/Documents…jeje)

import pandas as pd
import numpy as np
pd.set_option("display.max_rows", None, "display.max_columns", None)

# Lectura de datos
vinos = pd.read_csv('C:\\Users\\Guille\\Documents\\MineriaDatos_2022_23\\Datos\\DatosVino.csv')
vinos.head()

##    ID  Beneficio  Compra  Acidez  AcidoCitrico  Azucar  CloruroSodico  \
## 0   2        515       1    0.16         -0.81   26.10         -0.425   
## 1   4        585       1    2.64         -0.88   14.80          0.037   
## 2   8          0       0    0.29         -0.40   21.50          0.060   
## 3  11        775       1   -1.22          0.34    1.40          0.040   
## 4  12        596       1    0.27          1.05   11.25         -0.007   
## 
##    Densidad    pH  Sulfatos  Alcohol Etiqueta  CalifProductor Clasificacion  \
## 0   1.02792  3.38      0.70    144.0        M               2           ***   
## 1   0.99518  3.12      0.48     22.0        M               3           ***   
## 2   0.99572  3.49      1.21     10.3        R               3             ?   
## 3   1.03236  3.20       NaN     11.6        B               2           ***   
## 4   0.99620  4.93      0.26     15.0        R               1             ?   
## 
##    Region  PrecioBotella  
## 0     1.0           1.00  
## 1     3.0           3.38  
## 2     1.0           3.72  
## 3     2.0           6.23  
## 4     2.0           2.44

Este código está preparado para funcionar sobre los DatosVino. Se presenta en la siguiente tabla, la información sobre las variables contenidas en el archivo.

Fig1. Descripción de los datos En la tabla de información de variables que acompaña a los datos, podemos comprobar los tipos que deberían tener las variables y sus posibles limitaciones sobre rangos, categorías “oficiales” etc. Esto nos sirve de guía para ir chequeando lo que deberíamos tener y lo que realmente encontramos (dos cosas que a menudo no van de la mano precisamente…)

Que comprobar?

En el archivo total comprobaremos lo siguiente:

1- Tipos de variables. Todos los factores, lo son realmente en python?

2- Valores mal codificados. Hay missings no declarados tipo -1, 99999? Las categorías de las nominales son las que deben?

3- Valores fuera de rango. Hay alguna limitación sobre el rango de alguna variable que no se cumpla?

4- Variables nominales o Factores con categorías minoritarias. Las categorías con baja representación puede causar muchos problemas en los modelos por falta de base muestral para la estimación de los parámetros correspondientes a la pertenencia a esa categoría. Por ello, es conveniente echar un vistazo y recodificar las vairables uniendo categorías muy poco representativas con otras cuya unión tenga algún sentido (tienen comportamiento similar frente a la objetivo, la variable tiene caracter ordinal por lo que la unión con mayor sentido sería hacia categorías adyacentes..)

Con estas cosas ya arregladas, nos vamos a los dos grandes “caballos de batalla” de la depuración. Este proceso para la gestión de outliers y missings lo vamos a llevar a cabo sobre el archivo que llamaremos input y que contendrá solamente a los predictores sin incluir las variables objetivo continua y binaria. No es buena idea gestionar estas cosas en las objetivo puesto que las tomamos como la verdad verdadera y si tocamos algo ya estamos sesgando de alguna manera la variable de supervisión del proceso de modelización. Más adelante discutimos con mayor profundidad este hecho.

5- Outliers. Incidencia y tratamiento (pasar a missing, eliminar, winsorizar)

6- Missings. Incidencia y tratamiento (Eliminar, imputación simple por media, mediana, aleatorio, imputación por modelos)

Tipos de variables

Con la función info, podemos obtener un summary de los tipos de variables almacenados en el dataset como python lo ha entendido en la lectura.

# Información del dataset
vinos.info()

## <class 'pandas.core.frame.DataFrame'>
## RangeIndex: 6365 entries, 0 to 6364
## Data columns (total 16 columns):
##  #   Column          Non-Null Count  Dtype  
## ---  ------          --------------  -----  
##  0   ID              6365 non-null   int64  
##  1   Beneficio       6365 non-null   int64  
##  2   Compra          6365 non-null   int64  
##  3   Acidez          6365 non-null   float64
##  4   AcidoCitrico    6365 non-null   float64
##  5   Azucar          6365 non-null   float64
##  6   CloruroSodico   6066 non-null   float64
##  7   Densidad        6365 non-null   float64
##  8   pH              6170 non-null   float64
##  9   Sulfatos        5761 non-null   float64
##  10  Alcohol         6365 non-null   float64
##  11  Etiqueta        6365 non-null   object 
##  12  CalifProductor  6365 non-null   int64  
##  13  Clasificacion   6365 non-null   object 
##  14  Region          6258 non-null   float64
##  15  PrecioBotella   6365 non-null   float64
## dtypes: float64(10), int64(4), object(2)
## memory usage: 795.8+ KB

Se observa que no todas las variables tienen asignado el tipo correcto de datos. Identificamos factores como Compra (columna 3), Etiqueta (columna 12), Clasificación (columna 14) y Región (columna 15).

Valores únicos de las variables

Con nunique podemos obtener el número de valores distintos por cada variable para evaluar cuáles deberían ser factores.

# Número de valores distintos por variable
vinos.nunique()

## ID                6365
## Beneficio          983
## Compra               2
## Acidez             659
## AcidoCitrico       523
## Azucar            1639
## CloruroSodico     1406
## Densidad          3614
## pH                 436
## Sulfatos           543
## Alcohol            312
## Etiqueta            10
## CalifProductor      13
## Clasificacion        5
## Region               3
## PrecioBotella      598
## dtype: int64

En principio, claramente los factores seguros tienen menos de 10 valores distintos a excepción de Etiqueta, algo pasa ahí…Ya lo estudiaremos.

En este punto interesante charla sobre la posible dualidad continuo-categórica de la variable CalifProductor. Por una parte, tiene 13 valores, por lo que supera el umbral comentado de 10 valores distintos para ser considerada como tal. Por otra, con bajo número de valores, la linealidad con la respuesta se hace compleja, con lo que habría que comprobar la existencia de la misma. Considerando la variable como categórica, tenemos la ventaja de poder captar relaciones no lineales puesto que se modela la pertenencia a cada una de las categorías en relación a una de referencia. Sin embargo, hemos de ser conscientes del número de parámetros que consumirá en nuestro modelo \(k-1\) siendo K el número de niveles del factor.

Veamos su histograma.

import plotly.express as px
# Histograma de calificación del productor
fig = px.histogram(vinos, x="CalifProductor") 
fig.show()

Por tanto podemos adoptar la estrategia de mantener la variable continua por ese aspecto chi-cuadrado que no nos ha disgustado y, a su vez crear una variable categórica y posteriormente evaluar su distribución uniendo aquellas categorías minoritarias siempre con algún sentido, en este caso ordinal. Así tendremos las dos posibilidades para probar en los modelos de predicción que generemos.

Estas son posibles alternativas, en base a puro empirismo y un poco también porque se ve venir que dejar esa variable en un factor como pocas categorías que mantenga su poder predictivo es difícil tarea, decidimos mantenerla como continua.

Conversión a factor de toda variable con menos de 10 valores distintos

Utilizamos el .loc para filtrar por condición sobre valores únicos de variables y creamos la lista de columnas que deberían pasar a factor. Posteriormente, se utiliza astype para pasar a categórica.

# Lista de columnas con menos de 10 valores distintos. Potenciales factores!
to_factor = list(vinos.loc[:,vinos.nunique() < 10]);  

# Podemos cambiar el tipo de todas ellas a factor de una vez
vinos[to_factor] = vinos[to_factor].astype('category')

Vamos a contar el número de valores únicos de las variables numéricas por si nos hemos dejado algo por ahí.

# Número de valores distintos por variable
vinos.nunique()

## ID                6365
## Beneficio          983
## Compra               2
## Acidez             659
## AcidoCitrico       523
## Azucar            1639
## CloruroSodico     1406
## Densidad          3614
## pH                 436
## Sulfatos           543
## Alcohol            312
## Etiqueta            10
## CalifProductor      13
## Clasificacion        5
## Region               3
## PrecioBotella      598
## dtype: int64

En principio la única cosa rara es la variable Etiqueta que no debería tener tantos valores distintos..

Descirptivos para las variables

# Número de valores distintos por variable
vinos.describe()

##                  ID    Beneficio       Acidez  AcidoCitrico        Azucar  \
## count   6365.000000  6365.000000  6365.000000   6365.000000   6365.000000   
## mean    8010.702278   452.380204     0.331214      0.314350   4718.669780   
## std     4654.939139   308.380542     0.787534      0.861428  21192.546521   
## min        2.000000     0.000000    -2.790000     -3.240000   -127.100000   
## 25%     3980.000000   236.000000     0.130000      0.020000      0.900000   
## 50%     8065.000000   480.000000     0.280000      0.310000      5.000000   
## 75%    12027.000000   671.000000     0.650000      0.580000     22.600000   
## max    16128.000000  1568.000000     3.680000      3.860000  99999.000000   
## 
##        CloruroSodico     Densidad           pH     Sulfatos      Alcohol  \
## count    6066.000000  6365.000000  6170.000000  5761.000000  6365.000000   
## mean        0.051348     0.994204     3.202207     0.526659    16.250102   
## std         0.322715     0.026417     0.678330     0.948039    25.598217   
## min        -1.171000     0.888090     0.540000    -3.120000    -4.500000   
## 25%        -0.032750     0.988245     2.960000     0.280000     9.000000   
## 50%         0.046000     0.994400     3.190000     0.500000    10.500000   
## 75%         0.146750     1.000600     3.460000     0.880000    12.800000   
## max         1.351000     1.099240     6.050000     4.210000   150.000000   
## 
##        CalifProductor  PrecioBotella  
## count     6365.000000    6365.000000  
## mean         2.761508       2.610652  
## std          1.319127       1.480274  
## min          0.000000       1.000000  
## 25%          2.000000       1.420000  
## 50%          3.000000       2.190000  
## 75%          3.000000       3.440000  
## max         12.000000      11.440000

Distintas formas de echar el vistazo a las distribuciones de las variables, donde prestaremos atención al describe que nos informa sobre cuartiles y media, así como valores perdidos y máximos. Así, observamos que azucar tiene valores 99999 sospechosos y sulfatos 604 valores ausentes (NA), que alcohol debe tener distribución asimétrica positiva por valores posiblemente altos, de hecho es un % y no debería superar 100.

Podemos sacar los descriptivos básicos para las categóricas.

# Número de valores distintos por variable
vinos.describe(exclude=np.number)

##         Compra Etiqueta Clasificacion  Region
## count     6365     6365          6365  6258.0
## unique       2       10             5     3.0
## top          1        R            **     3.0
## freq      4998     2380          1754  2132.0

Inspección gráfica

Podemos inspeccionar las distribuciones gráficamente para completar la exploración. En este punto, nos podemos plantear la creación de un par de funciones gráficas que nos faciliten el trabajo, de tal forma que las podamos aplicar al dataset completo para visualizarlo de un plumazo.

En primer lugar, como me gusta ver el histograma y el boxplot en relación para cada variable, vamos a definir una función que nos saque ambos gráficos jusntos. Podríamos adoptar varias formas de pensar de cara a la programación, por flexibilidad, decidimos que la función trabaje sobre una potencial columna o variable, así podremos aplicarla a todas o parte de las columnas del dataset.

import seaborn as sns
import matplotlib.pyplot as plt

def histogram_boxplot(data, xlabel = None, title = None, font_scale=2, figsize=(9,8), bins = None):
    """ Boxplot and histogram combined
    data: 1-d data array
    xlabel: xlabel 
    title: title
    font_scale: the scale of the font (default 2)
    figsize: size of fig (default (9,8))
    bins: number of bins (default None / auto)

    example use: histogram_boxplot(np.random.rand(100), bins = 20, title="Fancy plot")
    """
    # Definir tamaño letra
    sns.set(font_scale=font_scale)
    # Crear ventana para los subgráficos
    f2, (ax_box2, ax_hist2) = plt.subplots(2, sharex=True, gridspec_kw={"height_ratios": (.15, .85)}, figsize=figsize)
    # Crear boxplot
    sns.boxplot(x=data, ax=ax_box2)
    # Crear histograma
    sns.histplot(x=data, ax=ax_hist2, bins=bins) if bins else sns.histplot(x=data, ax=ax_hist2)
    # Pintar una línea con la media
    ax_hist2.axvline(np.mean(data),color='g',linestyle='-')
    # Pintar una línea con la mediana
    ax_hist2.axvline(np.median(data),color='y',linestyle='--')
    # Asignar título y nombre de eje si tal
    if xlabel: ax_hist2.set(xlabel=xlabel)
    if title: ax_box2.set(title=title, xlabel="")
    # Mostrar gráfico
    plt.show()

histogram_boxplot(vinos.Beneficio, bins = 20, font_scale=1, title="Distribción de Beneficio")

Seguidamente vamos a implementar una función sencilla para sacar un barplot (aquí podéis indagar posibilidades de otros paquetes de visualización para customizar vuestros informes al gusto). Vamos a tirar de plotly a ver..


def cat_plot(col):
     if col.dtypes == 'category':
        fig = px.bar(col.value_counts())
        #fig = sns.countplot(x=col)
        return(fig)

# Aplicación a una variable en particular    
cat_plot(vinos.Compra)

cat_plot(vinos.Etiqueta)

Una vez definidas las funciones para gráficos de variables numéricas y nominales, podemos integrar el proceso en una sola función que pueda ser aplicada a todas o parte de las columnas y distinga el tipo de las mismas para escoger el gráfico adecuado y mostrarlo.


def plot(col):
     if col.dtypes != 'category':
        print('Cont')
        histogram_boxplot(col, xlabel = col.name, title = 'Distibución continua')
     else:
        print('Cat')
        cat_plot(col)
        


vinos.iloc[:,1:9].apply(plot)

## Cont
## Cat
## Cont
## Cont
## Cont
## Cont
## Cont
## Cont
## Beneficio        None
## Compra           None
## Acidez           None
## AcidoCitrico     None
## Azucar           None
## CloruroSodico    None
## Densidad         None
## pH               None
## dtype: object

Corrección de errores detectados

Etiqueta tiene 10 categorías y debería tener 5. Hacemos una tabla de frecuencias de las cetegorías para indagar.

# Tabla de frecuencias
vinos.Etiqueta.value_counts()

## R     2380
## M     1357
## B     1282
## r      420
## m      230
## MM     216
## b      209
## MB     191
## mm      40
## mb      40
## Name: Etiqueta, dtype: int64

Parece que hay mayúsculas y minúsculas. Habrá que unificar para que la variable tenga sentido. Lo más inmediato es pasar a mayúsculas o minúsculas todo y sencillamente se arreglará.

Podemos utilizar el apply con lambda para ver cómo qudaría el cambio (también podríamos utilizarlo para cambiarlo asignando a la variable y “pisando” la información)

# Utilizar apply con lambda para ver como quedaría el cambio
vinos['Etiqueta'].apply(lambda x: x.upper()).value_counts(sort=False)

## M     1587
## R     2800
## B     1491
## MB     231
## MM     256
## Name: Etiqueta, dtype: int64

Tiene buena pinta. Podemos hacer el cambio sencillamente con operaciones sobre cadenas de caracteres de pandas .str y la función upper() para transformación en mayúsculas.

# Realizar el cambio en la propia variable y cambiar el tipo a categórica
vinos['Etiqueta'] = vinos['Etiqueta'].str.upper().astype('category')

Ya que estamos con motivación, vamos a reaordenar el factor para que MM sea el primer nivel o el más bajo y MB sea el último, manteniendo la naturaleza ordinal de la variable.

# Vamos a ordenar el factor como nos gustaría tenerlo
vinos["Etiqueta"] = vinos["Etiqueta"].cat.reorder_categories(['MM','M','R','B','MB'])

# A ver como está el tema ahora
vinos.Etiqueta.value_counts()

## R     2800
## M     1587
## B     1491
## MM     256
## MB     231
## Name: Etiqueta, dtype: int64

La tabla de frecuencia de value_counts() ordena por defecto las categorías por volumen de registros o frecuencia, pero siempre podemos decirle que no lo haga, con lo que mostrará el órden original del factor.

# Tabla de frecuencias sin ordenar por volumen
vinos.Etiqueta.value_counts(sort=False)

## MM     256
## M     1587
## R     2800
## B     1491
## MB     231
## Name: Etiqueta, dtype: int64

Azucar. Tiene un gráfico directamente ilegible…está claro que hay un NA no declarado con valor 99999. Hay que reemplazarlo por un verdadero missing. En python, podemos asignar un valor missing con la función nan de numpy. Al realizar el cambio con la filosofía inplace, se pisa la información original en la columna del dataset.

Cuidado con utilizar esto cuando no tenemos muy claro que el cambio es correcto!! Habría que volver atrás a datos brutos y volver a ejecutar.. que bueno, por otro lado es gratis, pero mejor para pruebas trabajar con copias o con inplace=False.

# Quitamos ese 99999 de la variable Azucar y lo pasamos a NA
vinos.Azucar.replace(99999,np.nan,inplace=True)

# Comprobamos el nuevo máximo
vinos.Azucar.max()

## 141.15

Alcohol. La tabla de info de las variables nos dice que es un % y debe etar por tanto restingido al rango 0-100. Todo valor fuera de ese rango habrá que tomarlo como valor perdido pues no podemos asignar claramente a otros valores en ppio.

# Corregimos valores fuera de rango de Alcohol (recordemos que es un %). Opción con np.where de nuumpy
#vinos.Alcohol = np.where((vinos['Alcohol'] < 0) | (vinos['Alcohol'] >100), np.nan,vinos.Alcohol)

# Opción con .loc de pandas
vinos.loc[~vinos.Alcohol.between(0, 100), "Alcohol"] = np.nan

# Comprobamos nuevos límites de la variable
vinos.Alcohol.min(),vinos.Alcohol.max()

## (0.0, 26.5)

Clasificación. Veíamos esa categoría ? que nos llama la atención. Veamos si se trata de algo residual y podemos pasar a missing para imputar posteriormente o si es una categoría con importancia en volumen y patrón frente al objetivo.

# Frecuencias de las categorías de clasificación
vinos.Clasificacion.value_counts()

# Tabla de contingencia con la variable objetivo binaria Compra

## **      1754
## ?       1680
## *       1535
## ***     1074
## ****     322
## Name: Clasificacion, dtype: int64

pd.crosstab(index=vinos['Compra'], columns=vinos['Clasificacion'])

## Clasificacion     *    **   ***  ****     ?
## Compra                                     
## 0               306    46     0     0  1015
## 1              1229  1708  1074   322   665

Reemplazar los ? por Desconocido pues representa una categoría importante en volumen y con un posible patrón interesante ante la variable objetivo. Se intuye que los vinos de clasificación desconocida tienen bastante menor probabilidad de compra!

# Reemplazar los ? por Deconocido pues representa una categoría importante en volumen y con un posible patrón interesante
vinos.Clasificacion.replace('?','Desc',inplace=True)

Una vez libres de errores graves, las variables están preparadas para la gestión de outliers y missing. Para ello, es importante separa las variables objetivo y trabajar en el archivo de predictores. No es habitual tocar las variables objetivo puesto que representan nuestra verdad verdadera, son las variables de supervisión y se presuponen bien recogidas.

Imaginemos que se nos presenta la situación en la que tenemos valores missing en las objetivo, que deberíamos hacer? Pues tratar estas instancias como un conjunto de test sobre el que podríamos hacer predicciones y valorar si el modelo parece tener sentido. El problema es que no podremos evaluar la calidad de las estimaciones mediante el error cometido puesto que no tenemos su verdad verdadera.

#Indico la variableObj, el ID y las Input 
# los atípicos y los missings se gestionan sólo de las input
varObjCont = vinos.Beneficio
varObjBin = vinos.Compra
imput = vinos.drop(['Beneficio','Compra'],axis=1)

imput.head()

##    ID  Acidez  AcidoCitrico  Azucar  CloruroSodico  Densidad    pH  Sulfatos  \
## 0   2    0.16         -0.81   26.10         -0.425   1.02792  3.38      0.70   
## 1   4    2.64         -0.88   14.80          0.037   0.99518  3.12      0.48   
## 2   8    0.29         -0.40   21.50          0.060   0.99572  3.49      1.21   
## 3  11   -1.22          0.34    1.40          0.040   1.03236  3.20       NaN   
## 4  12    0.27          1.05   11.25         -0.007   0.99620  4.93      0.26   
## 
##    Alcohol Etiqueta  CalifProductor Clasificacion Region  PrecioBotella  
## 0      NaN        M               2           ***    1.0           1.00  
## 1     22.0        M               3           ***    3.0           3.38  
## 2     10.3        R               3          Desc    1.0           3.72  
## 3     11.6        B               2           ***    2.0           6.23  
## 4     15.0        R               1          Desc    2.0           2.44

Valores atípicos

Para facilitarnos la vida y complementar la idea que tenemos ya sobre las distribuciones de las variables, llevamos a cabo un conteo de los valores que se consideran extremos según un consenso de dos criterios distintos. En primer lugar, se distingue variable simétrica o posiblemente no, para aplicar media + 3 sd ó mediana + 8 mad, respectivamente. Recordamos en este punto que todas las medidas de dispersión basadas en la mediana o cuartiles son muy poco sensibles a la presencia de asimetría en la distribución, siendo por ello más fiables en este caso. Por otro lado, aplicamos el clásico criterio del boxplot umbrales en cuartil1 - 3IQR y cuartil3+ 3IQR.

Para aclarar, veamos un momento las asimetrías de las variables numéricas.

vinos.select_dtypes(include=np.number).apply(lambda x: x.skew())

## ID                0.003339
## Beneficio         0.023872
## Acidez            0.029079
## AcidoCitrico     -0.024261
## Azucar            0.022391
## CloruroSodico     0.012041
## Densidad         -0.009353
## pH                0.003311
## Sulfatos         -0.062373
## Alcohol           0.267099
## CalifProductor    1.663553
## PrecioBotella     1.120795
## dtype: float64

Asimetrías en valor absoluto mayores a la unidad son signo de distribución significativamente sesgada a la derecha/positiva (+) o izquierda/negativa (-).

Antes de gestionar aquellos valores detectados como outliers, valoramos la incidencia en cada variable contando el porcentaje de registros atípicos en relación al total de registros. Con ello, evaluamos si estamos ante un problema grave (tener más de un 20% de valores atípicos en alguna variable es síntoma de distribución bimodal, hay en realidad dos poblaciones jugando en la distribución) o si por el contrario, se trata de un % bajo con el que podamos lidiar sin graves consecuencias en las distribuciones de las variables.

Vamos a crear una función gestiona_outliers para facilitar el trabajo en este sentido. La idea es tener una sola función que implemente:

1- Visualización de la incidencia de outliers –> Toma de decisiones

2- Gestión de los mismos:

2.1- Winsorizar: Colapsar todo valor mayor o menor que los umbrales de cierto rango (dado por percentiles habitualmente) para igualarlo con esos umbrales. Como ya comentábamos, este método conserva la integridad del datos en el sentido de que si era muy grande por la derecha ahora estará en el extremo, pero puede provocar carga excesiva de las colas de la distribución, generando quizá algunas rarezas en picos de densidad..

2.2- Convertir en NAs: Asumiendo que los valores extremos son inaceptables (es decir, se sospecha que no puede ser tan alto o bajo por alguna razón). Podríamos pensar en que sea un fallo y por tanto decidir pasar a missing esos valores para luego gestionarlos con los valores perdidos. Esta técnica podría incurrir en errores graves de cambio del sentido de la información. Por ejemplo, en este caso un vino muy caro, pero que es un vino real y ese es su precio, es lo que hay…bueno considerándo que es erróneo, pasamos a Na e imputamos por media o valor aleatorio, quedando este valor de precio en un entorno central…tal vez hay características que generaban un patrón de precio alto que ahora serán mal aprendidas por el modelo… Hay que andarse con ojo con estas cosas..

En cualquier caso, ahí queda la función para utilizarla como se considere y tal vez incluir o desarrollar nuevas formas de tratamiento..


from scipy.stats.mstats import winsorize

# Función para gestionar outliers
def gestiona_outliers(col,clas = 'check'):
     print(col.name)
     # Condición de asimetría y aplicación de criterio 1 según el caso
     if abs(col.skew()) < 1:
        criterio1 = abs((col-col.mean())/col.std())>3
     else:
        criterio1 = abs((col-col.median())/col.mad())>8
     
     # Calcular primer cuartil     
     q1 = col.quantile(0.25)  
     # Calcular tercer cuartil  
     q3 = col.quantile(0.75)
     # Calculo de IQR
     IQR=q3-q1
     # Calcular criterio 2 (general para cualquier asimetría)
     criterio2 = (col<(q1 - 3*IQR))|(col>(q3 + 3*IQR))

     lower = col[criterio1&criterio2&(col<q1)].count()/col.count()
     upper = col[criterio1&criterio2&(col>q3)].count()/col.count()
     
     # Salida según el tipo deseado
     if clas == 'check':
            return(lower*100,upper*100,(lower+upper)*100)
     elif clas == 'winsor':
            return(winsorize(col,(lower,upper)))
     elif clas == 'miss':
            print('\n MissingAntes: ' + str(col.isna().sum()))
            col.loc[criterio1&criterio2] = np.nan
            print('MissingDespues: ' + str(col.isna().sum()) +'\n')
            return(col)
          
# Llamada en modo check
imput.select_dtypes(include=np.number).copy().apply(lambda x: gestiona_outliers(x))

## ID
## Acidez
## AcidoCitrico
## Azucar
## CloruroSodico
## Densidad
## pH
## Sulfatos
## Alcohol
## CalifProductor
## PrecioBotella
##     ID    Acidez  AcidoCitrico    Azucar  CloruroSodico  Densidad        pH  \
## 0  0.0  0.816968      0.848390  0.956307       0.807781  0.942655  0.794165   
## 1  0.0  0.879811      0.785546  0.972795       0.873722  0.942655  0.696921   
## 2  0.0  1.696779      1.633936  1.929101       1.681503  1.885310  1.491086   
## 
##    Sulfatos   Alcohol  CalifProductor  PrecioBotella  
## 0  1.024128  0.000000        0.000000            0.0  
## 1  0.815831  0.317885        0.094266            0.0  
## 2  1.839958  0.317885        0.094266            0.0

Como pequeño resumen, lo que buscamos con el tratamiento de outliers es identificar valores verdaderamente extremos y que no sea típicos, lo cual necesariamente implica que han de ser pocos!! Y esto es muy relevante. Si son el 20% de los registros…pues tal vez estamos ante dos poblaciones distintas… habría que abordar el problema de otra forma..tal vez tratar cada población por separado.

# Crear copia para evitar pisar información
vinCont = imput.select_dtypes(include=np.number).copy()

# Aplicar la gestión de outliers en modelo winsor
vinos_winsor = vinCont.apply(lambda x: gestiona_outliers(x,clas='winsor'))

# Contemos si ha desaparecido algún resgitro o algo 
# vinos_winsor.apply(lambda x: x.isna().sum()/x.count()*100)

# Veamos los valores mínimos de acidez ahora
# vinos_wins.sort_values(by='Acidez').head()

# Veamos los valores mínimos de acidez ahora
# vinos.sort_values(by='Acidez').head()

## ID
## Acidez
## AcidoCitrico
## Azucar
## CloruroSodico
## Densidad
## pH
## Sulfatos
## Alcohol
## CalifProductor
## PrecioBotella

Recordemos que la búsqueda de outliers atañe en exclusiva a las variables numéricas por lo que el archivo de salida no tiene los factores. Vamos a juntar numéricas y factores en un solo archivo para valorar la incidencia de perdididos en el conjunto completo.


# Juntar columnas con join
imput_wins = vinos_winsor.join(vinos.select_dtypes(exclude=np.number))

Muchas alternativas disponibles, lo importante es conocer las fortalezas y debilidades de cada uno y aplicarlo con lógica según resulte más conveniente para los datos que manejamos.

Valores perdidos

Entramos de lleno en la segunda gran gestión que debemos llevar a cabo antes de la modelización. Los archiconocidos valores perdidos. En primer lugar comentar que, llegados a este punto tenemos valores perdidos de dos fuentes distintas, por una parte los que vienen “de serie” en el dataset que podemos asociar a falta de medida en la recogida de los datos y aquí existe toda una teoría sobre los mecanismos de aparición de dichos missing, completamente al azar (MCAR), al azar(MAR) o nada de azar y existe patrón (MNAR).

Dejo por aquí un poco más de información y diversos métodos propuestos para la imputación de estos valores.

https://stefvanbuuren.name/fimd/sec-MCAR.html

Nuestro objetivo, en este limitado módulo, será valorar la incidencia de los valores perdidos y conocer los métodos usuales de imputación univariante (cada variable independientemente de las otras) con sus pros y contras.

Que podemos hacer con los missings?

Eliminar u omitir del dataset

Eliminar “por lista” los Nas, es decir, eliminar del dataset toda observación con al menos un NA en al menos 1 variable. Se consigue liberar el dataset de perdidos (algunas técnicas y análisis no permiten su inclusión) pero se puede incurrir en los problemas que comentamos en el siguiente punto.

Nada

Podemos obviar la presencia de valores perdidos y ya el modelo se encargará de quitarlos “por lista”, es decir, eliminar del análisis toda observación con al menos un NA en el set de variable escogido (importante este punto! si elegimos un modelo con las variables v1:v5, no se eliminarán los registros completos en estas variables pero con perdidos en v6:v10, he aquí la diferencia con la omisión de missings de primeras para toto el archivo). Esto es habitual y se puede hacer, ahora bien no está exento de peligros. Veamos.

Peligro de los missings cruzados. Imaginamos el caso en que tenemos 10 variables y 100 registros y cada una de ellas tiene un 5% de perdidos…No parece mucho con lo cual los quitamos por lista. Nuestro pensamiento es 5% pero, si se da el caso de que los NA de cada variable aparecen en registros distintos…entonces tenemos 5*10= 50% de los registros con al menos un perdido…hemos perdido la mitad de la información!!!
Sesgo por valores perdidos. El simple hecho de eliminar observaciones por el mero hecho de que presenten perdidos puede introducir un importante sesgo de selección de registros en los modelos. Imaginemos que la gente mayo tiende a no contestar ciertas preguntas en una encuesta, eliminamos y nos quedamos con muy poca gente mayor por lo que las conclusiones con toda seguridad estarán sesgadas hacia los jóvenes.

Imputar sus valores.

No queremos exponernos a lo anterior por lo que se puede adoptar la estrategia de asignar un valor a estos datos no conocidos.

Imputación por valores centrales (media, mediana): Muy habitual asignar valores centrales ya que no alteran la distribución de las variables como tal. El gran inconveniente de este método es a subestimación de la verdadera varianza de la variable ya que estamos centrando demasiado la distribución haciendo que artificialmente su varianza se reduzca, en ocasiones muy drásticamente.
Imputación por valores aleatorios: Si no queremos centrar tanto la distribución, podemos optar por asignar al azar valores observados de las distribuciones de cada variable a los registros con NA. De esta forma, cualquier valor en el rango observado puede ser asignado a los faltantes. El gran inconveniente de esto es la dependencia del azar.
Imputación por modelos multivariantes: Muchas opciones en este apartado. Existen métodos que tienen en cuenta los valores observados de otras variables para asignar el valor más “plausible” a la variable perdida en un sentido conjunto. Una de las alternativas es generar imputaciones por un modelo de regresión por ejemplo, así para imputar Azucar utilizaremos un modelo que estime los valores de azucar en base a las demás variables (que se ajustará con los valores validos por lista) y posteriormente predeciremos los perdidos de azucar mediante este modelo generado. De esta misma forma existen modelos de imputación por random forest (missforet), vecinos más cercanos (knn), cadenas de Markov (hmisc, amelia) y gran cantidad de aproximaciones de imputación múltiple (imputar n veces y promediar). El mayor problema de estos métodos suele radicar en el posible sobreajuste a los datos de training.

El primer paso es valorar la incidencia de los valores perdidos ya que si no representan gran proporción, no existe gran peligro de cambio en la distribución de las variables con independencia del método utilizado para la imputación.

Incidencia

Vamos a actuar por dos vías, filas y columnas para contar el % de perdidos en cada registro y en cada variable. Con esta información tomaremos decisiones sobre gestión.

Missings por variable.

#Proporción de missings por variable 
imput_wins.apply(lambda x: x.isna().sum()/x.count()*100)

## ID                 0.000000
## Acidez             0.000000
## AcidoCitrico       0.000000
## Azucar             4.946414
## CloruroSodico      4.929113
## Densidad           0.000000
## pH                 3.160454
## Sulfatos          10.484291
## Alcohol            6.491551
## CalifProductor     0.000000
## PrecioBotella      0.000000
## Compra             0.000000
## Etiqueta           0.000000
## Clasificacion      0.000000
## Region             1.709811
## dtype: float64

Controlado esto, vemos que sulfatos es la variable más peligrosa con más de un 10% de perdidos..Digamos que es la variable en la que mayores dudas tenemos respecto a la conservación de su distribución tras la imputación. Por lo demás, valores relativamente bajos.

Missings por observación.

Calcularemos ahora el % de missings por observación y vamos a aplicar un truquiflor que a veces nos ayuda mucho a controlar las imputaciones. Se trata de generar una variable en el archivo que cuanta la proporción de perdidos que tiene cada registro. De esta forma siempre tenemos una huella de los registros con alta carga de imputaciones por si necesitamos saber algo sobre ellos en la etapa de modelado. Es más, esta será una variable que incluiremos en el modelo para valorar si puede generar un patrón de comportamiento respecto a la variable objetivo.

#Proporción de missings por observación (como una nueva columna del dataset)
imput_wins['prop_missings'] = imput_wins.apply(lambda x: x.isna().sum()/x.count()*100,axis=1)

# Valoramos distribución
imput_wins.prop_missings.describe()

## count    6365.000000
## mean        2.165129
## std         3.944149
## min         0.000000
## 25%         0.000000
## 50%         0.000000
## 75%         7.142857
## max        25.000000
## Name: prop_missings, dtype: float64

Vamos a ordenar el archivo por la nueva variable creada para ver el aspecto.

imput_wins.sort_values(by='prop_missings', ascending=False).head()

##          ID  Acidez  AcidoCitrico  Azucar  CloruroSodico  Densidad    pH  \
## 2781   6990    1.40          0.30     NaN            NaN   1.04263  4.00   
## 1344   3341    0.26          0.27    18.2          0.048   1.04227   NaN   
## 5402  13598   -0.45          0.26     NaN            NaN   0.99907  2.21   
## 4385  11020    0.27         -0.01     NaN          0.295   0.99566   NaN   
## 1352   3357    1.02          1.51     NaN            NaN   0.99300  3.37   
## 
##       Sulfatos  Alcohol  CalifProductor  PrecioBotella Compra Etiqueta  \
## 2781       NaN     15.1               3           1.30      1        R   
## 1344       NaN      3.8               1           2.62      1        B   
## 5402      0.68      NaN               3           3.10      1        M   
## 4385     -0.46      NaN               3           2.35      1        R   
## 1352      0.55     11.5               2           1.60      0       MB   
## 
##      Clasificacion Region  prop_missings  
## 2781            **    2.0           25.0  
## 1344            **    NaN           25.0  
## 5402          Desc    1.0           25.0  
## 4385            **    2.0           25.0  
## 1352          Desc    NaN           25.0

El siguiente código pretende eliminar registros y observaciones con más de la mitad de su información perdida puesto que resulta arriesgado imputar tal cantidad de datos perdidos. Es evidente que en nuestro caso no aplica y si lo ejecutamos no habrá cambio alguno.

Es importante saber que tenemos que aplicar el mismo filtro a las variables objetivo para que al unir el input depurado con ellas, los registros cuadren!

#elimino las observaciones y las variables con más de la mitad 
# de datos missings (si no hay ninguna, no ejecuto este código)
input %>% filter(prop_missings< 0.5) %>% select(!(names(
  prop_missingsVars)[prop_missingsVars>0.5]))-> imput 

#Actualizar las observaciones de las variables objetivo
varObjBin<-varObjBin[input$prop_missings<0.5] 
varObjCont<-varObjCont[input$prop_missings<0.5]

Coexistencia y patrones de missings

Podemos echar un vistazo a la correlación en la existencia de missings para valorar si existe algún patrón de aparición de los mismos. En caso de que observemos patrones de aparición, podemos tirar del hilo e indagar el porqué de ese comportamiento para decidir el método más adecuado para imputar.

Hay un paquete disponible llamado missingno que implementa un par de funciones gráficas interesantes para la detección de patrones de coaparición de los valores perdidios en diferentes variables.

#conda config --add channels conda-forge
#conda install missingno
# pip install missingno
import missingno as msno   

# Plot correlation heatmap of missingness
msno.matrix(imput_wins.sort_values(by='Azucar'))

No se intuyen patrones de aparición de missings con facilidad.

Veamos esto en el mapa de calor.

msno.heatmap(imput_wins)

Relaciones muy débiles con ese máximo de 0.1 para pH con Región. Nada destacable.

Imputaciones

Llegados a este punto, hemos comprobado que la incidencia de valores perdidos no es preocupante en general (con sulfatos con ese 10% que nos invita a estar atentos) por lo que vamos a decidirnos por la imputación de los datos de cara a mantener la mayor base muestral posible para el ajuste de los futuros modelos a los datos.

Nos puede la curiosidad y nos planteamos que pasaría si elimináramos los NAs por lista del dataset. Cuantas observaciones perderíamos? Las variables que mayor carga de perdidos presentan, son realmente relevantes para un posible modelo?

imput_wins.dropna().describe()

##                  ID       Acidez  AcidoCitrico       Azucar  CloruroSodico  \
## count   4692.000000  4692.000000   4692.000000  4692.000000    4692.000000   
## mean    7996.172421     0.328282      0.324333     5.938544       0.055412   
## std     4671.828755     0.766658      0.832756    34.233343       0.321277   
## min        4.000000    -2.050000     -2.240000   -96.600000      -0.910000   
## 25%     3917.750000     0.130000      0.030000    -1.700000      -0.029000   
## 50%     8071.000000     0.280000      0.310000     4.100000       0.046000   
## 75%    12026.250000     0.650000      0.580000    16.125000       0.152250   
## max    16111.000000     2.680000      2.880000   140.650000       1.351000   
## 
##           Densidad           pH     Sulfatos      Alcohol  CalifProductor  \
## count  4692.000000  4692.000000  4692.000000  4692.000000     4692.000000   
## mean      0.994295     3.196036     0.533186    10.586679        2.762788   
## std       0.025540     0.668806     0.932729     3.498461        1.309899   
## min       0.915440     1.230000    -2.150000     0.000000        0.000000   
## 25%       0.988540     2.960000     0.290000     9.000000        2.000000   
## 50%       0.994400     3.190000     0.500000    10.400000        3.000000   
## 75%       1.000705     3.460000     0.880000    12.400000        3.000000   
## max       1.072380     6.050000     4.110000    26.500000       10.000000   
## 
##        PrecioBotella  prop_missings  
## count    4692.000000         4692.0  
## mean        2.610217            0.0  
## std         1.482378            0.0  
## min         1.000000            0.0  
## 25%         1.420000            0.0  
## 50%         2.180000            0.0  
## 75%         3.430000            0.0  
## max        11.440000            0.0

Nos quedamos con 4692 registros en total para todas las variables, lo que supone un 73% de la información del archivo. Nos hemos cargado el 27% de registros..

Como veremos cuando estudiemos las relaciones con las variables objetivo, en este dataset, las variables numéricas o características químicas (todas esas distribuciones super apuntadas) no generan especial patrón con las objetivo por lo que no serán predictores relevantes. Si nos cargamos ahora información en base a sus faltantes estaremos limitando o sesgando la capacidad del modelo (que seguramente no contenga a estas variables y pueda utilizar esos registros) por disminución (además no aleatoria) de base muestral.

Conclusión, imputemos!

Vamos a explorar algunas posibilidades de los paquetes sklearn y feature_engine para imputaciones simples y multivariantes de varios tipos.

Definiremos los imputadores disponibles y posteriormente los aplicaremos a los datos.


import sklearn.impute as skl_imp
from sklearn.experimental import enable_iterative_imputer
# Moda: Solo nominales
imputer_moda = skl_imp.SimpleImputer(strategy='most_frequent', missing_values=np.nan)
# knn: Solo numéricas
imputer_knn = skl_imp.KNNImputer(n_neighbors=3)
# Chain equations: solo numéricas
imputer_itImp = skl_imp.IterativeImputer(max_iter=10, random_state=0)

import feature_engine.imputation as fe_imp
# Aleatoria: numéricas y nominales
imputer_rand = fe_imp.RandomSampleImputer()
# Mediana: solo numéricas
imputer_median = fe_imp.MeanMedianImputer(imputation_method='median')
# Media: solo nominales
imputer_mean = fe_imp.MeanMedianImputer(imputation_method='mean')

Separaremos el dataset de nuevo en continuas y categóricas para aplicar los métodos que correspondan.

imput_wins_cont = imput_wins.select_dtypes(include=np.number)
imput_wins_cat = imput_wins.select_dtypes(exclude=np.number)

Posibilidades para las numéricas:

1- Nivel univariante

Comenzamos proponiendo métodos de imputación univariante que se basan en la información exclusiva de la distribución de la propia variable sin mirar a sus compañeras. En este plano, opciones ya comentadas como media, mediana o aleatorio.

# Media
vinos_winsor_mean_imputed = imputer_mean.fit(imput_wins_cont).transform(imput_wins_cont)
# Mediana
vinos_winsor_median_imputed = imputer_median.fit(imput_wins_cont).transform(imput_wins_cont)

2- Nivel multivariante

Podemos explorar un par de posibilidades que se utilizan en el mundillo. Por una parte, las imputaciones basadas en la técnica de los k vecinos más cercanos (digamos que es un método espacial) que asocia valores promedio de los k vecinos más cercanos a la observación perdida pero jugando en el espación R^k siendo k el número de variables. Es decir, identifica registros parecidos en general en todas las variables (que estén cerca en ese hiperespacio), lo cual es una buena idea. El mayor problema que tiene es la dependencia del valor k (no hay consenso en este aspecto y depende de los datos) y el posible sobreajuste a los datos de training como cualquier otro método multi.

# Fit/transform
imput_wins_knn_imputed = pd.DataFrame(imputer_knn.fit_transform(imput_wins_cont),columns=imput_wins_cont.columns)
imput_wins_itImp_imputed = pd.DataFrame(imputer_itImp.fit_transform(imput_wins_cont),columns=imput_wins_cont.columns)

Posibilidades para las nominales:

Para las variables categóricas podemos utilizar moda (la categoría más representada) o aleatorio.

imput_wins_moda_imputed = pd.DataFrame(imputer_moda.fit_transform(imput_wins_cat),columns=imput_wins_cat.columns)

Aplicar random a nivel general ya que acepta todo tipo.

imput_wins_rand_imputed = imputer_rand.fit(imput_wins).transform(imput_wins)
imput_wins_rand_imputed.describe()

##                  ID       Acidez  AcidoCitrico       Azucar  CloruroSodico  \
## count   6365.000000  6365.000000   6365.000000  6365.000000    6365.000000   
## mean    8010.702278     0.330902      0.315136     5.801131       0.051932   
## std     4654.939139     0.768507      0.841088    33.980154       0.318607   
## min        2.000000    -2.050000     -2.240000   -96.600000      -0.910000   
## 25%     3980.000000     0.130000      0.020000    -1.600000      -0.032000   
## 50%     8065.000000     0.280000      0.310000     4.000000       0.046000   
## 75%    12027.000000     0.650000      0.580000    15.900000       0.147000   
## max    16128.000000     2.680000      2.880000   141.150000       1.351000   
## 
##           Densidad           pH     Sulfatos      Alcohol  CalifProductor  \
## count  6365.000000  6365.000000  6365.000000  6365.000000     6365.000000   
## mean      0.994207     3.204671     0.542364    10.625656        2.760094   
## std       0.025614     0.671902     0.932792     3.547211        1.310444   
## min       0.915440     1.230000    -2.150000     0.000000        0.000000   
## 25%       0.988245     2.960000     0.290000     9.000000        2.000000   
## 50%       0.994400     3.190000     0.500000    10.400000        3.000000   
## 75%       1.000600     3.470000     0.890000    12.400000        3.000000   
## max       1.072380     6.050000     4.210000    26.500000       10.000000   
## 
##        PrecioBotella  prop_missings  
## count    6365.000000    6365.000000  
## mean        2.610652       2.165129  
## std         1.480274       3.944149  
## min         1.000000       0.000000  
## 25%         1.420000       0.000000  
## 50%         2.190000       0.000000  
## 75%         3.440000       7.142857  
## max        11.440000      25.000000

Toma de decisiones y guardado del archivo depurado

# Agregar variables objetivo al input ya limpio
vinos_wins_rand_imputed = pd.concat([imput_wins_rand_imputed, varObjCont,varObjBin], axis=1)

# Guardar archivo
vinos_wins_rand_imputed.to_csv('C:\\Users\\Guille\\Documents\\MineriaDatos_2022_23\\PARTE I_Depuracion y Regresiones\\Dia1_MDDepuracion\\DatosVinoDep_winsRand.csv')

Ya tenemos los datos depurados para poder empezar con el modelado. Es importante saber que a la hora de modelar utilizaremos este nuevo conjunto datosVinoDep_metodos y no el original, que para eso nos lo hemos trabajado.

Notas:

Debido a la elección del método de imputación por valores aleatorios, cada ejecución de este script generará un conjunto de datos depurados con ligeras diferencia en distribución de las variables sometidas a dicho proceso. Esto es natural y, dado que la verdadera distribución de los valores de los perdidos se deconoce en teoría, cualquiera de las soluciones sería plausible. Este comportamiento provoca que el posterior proceso de modelización con base en los datos depurados presente valores dependientes de la ejecución concreta de este primer paso. La idea es que nuestra limpieza de los datos sea relativamente robusta, por lo que los cambios deberían ser muy ligeros y no producirse cosas como cambio radical de influencia de las variables en los modelos o muy distintas capacidades de predicción sobre los datos de test. Si esto pasa, sospecharemos que nuestro proceso preliminar produce sesgos peligrosos y tendremos que repasar en especial el tratamiento de outliers y su conversión a valores perdidos (son lo que más influencia tienen en este hecho!)
Es evidente que muchas tipos de tratamiento de variables son posibles y normalmente es buena política “jugar” lo máximo con ellos y proponer distintas recodificaciones para aumentar nuestras posibilidades en la búsqueda de esos patrones ocultos que relacionan la variabilidad de las variables y que serán descubiertos por el modelo para ofrecer predicciones decentes.
Antes de lanzarse a las imputaciones es conveniente pensar si existen relaciones entre variables de manera natural, por ejemplo cuando dos continuas tienen una relación lineal fortísima digamos 0.9 de coeficiente de correlación r de Pearson, podría ser muy buena idea asignar directamente valores faltantes de una de ellas cuando la otra presenta valor conocido, así podemos reducir la incertidumbre y el sesgo de las imputaciones. De igual forma, relaciones transitivas tipo A ~ B*C (el producto de B y C presenta una cotrrelación alta con la variable A). Oye pues cuando A está perdida y ese registro tiene los valores de B y C conocido, en lugar de imputarle su valor de A, asignamos directamente.
En cuanto a código y ejecuciones en RMD, por las dudas. Los chunks contienen el código y pueden ser ejecutados por separado en su totalidad con el boton del play (parte superior derecha del chunk), línea a línea con ctrl+enter o el botón Run. Hemos de saber que a la hora de compilar el documento, las rutas deben ser válidas y contener los elementos que se requieren! El tipo de documento que se genera se controla en la cabecera y hay muchísimas posibilidades tanto en html como en word o pdf. Os animo a investigar formatos de vuestro agrado.